Nonthermal symmetry broken states in the strongly interacting Hubbard model 
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We study the time evolution of the antiferromagnetic order parameter after interaction quenches 
in the Hubbard model. Using the nonequilibrium dynamical mean field formalism, we show that 
the system, after a quench from intermediate to strong interaction, is trapped in a nonthermal state 
which is reminiscent of a photo-doped state and protected by the slow decay of doublons. If the 
effective doping of this state is low enough, it exhibits robust antiferromagnetic order, even if the 
system is highly excited and the thermal state thus expected to be paramagnetic. We comment on 
the implication of our findings for the stability of nonthermal superconducting states. 
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The rapid progress in the development of time-resolved 
spectroscopies enables direct measurements of dynamical 
symmetry breaking or the melting of long-range order 
in correlated materials. 1 Ultra- fast demagnetization af- 
ter femto-second optical pulses was first observed in fer- 
romagnetic Nii2 The transient closing of the charge den- 
sity wave (CDW) gap and excitations of CDW amplitude 
oscillations in lT-TaS 2 3 ^ and TbTe 3 ? were extensively 
studied using time-resolved photoemission. These exper- 
iments showed that the CDW gap disappears within less 
than 100 femto-seconds (essentially the time-resolution of 
the experiment). Ref. 9 reported the melting of orbital 
order and related structural transitions in Manganites af- 
ter a strong laser pulse. While the above studies focus 
on the destruction of some long-range ordered state (or 
at least the associated gap in the electronic spectrum) 
by a strong excitation, and emphasize the fast time-scale 
on which this process happens, there have also been re- 
cent experiments on photo-stimulated cuprates, which in- 
dicate the emergence of a nonthermal superconducting 
state which appears to be remarkably stabled 

Since all of these experiments involve transition metal 
compounds, they raise fundamental questions about the 
dynamics of symmetry breaking transitions and the sta- 
bility of symmetry broken states in strongly interacting 
electron systems. While the concept of a phase tran- 
sition is well understood in equilibrium, it is not obvi- 
ous how to apply these ideas if a transition occurs un- 
der nonequilibrium conditions. For example, it is un- 
clear whether and how the timescale for the melting of 
(quasi) long-range order is related to equilibrium corre- 
lation times. A recent theoretical study of a dynam- 
ical Kosterlitz-Thouless transition found an evolution 
through nonthermal (superheated) states after a quench, 
and showed that the transition to a disordered state can 
become very slow if the latter is close to the equilibrium 
phase transition.— This example indicates that even a 
strongly excited nonequilibrium state can differ from a 
high-temperature disordered state on long timescales. 



In this paper, we explore the dynamics of a symmetry 
broken state in the most fundamental model for corre- 
lated electron materials, the Hubbard model, after an in- 
teraction quench in the strongly correlated regime. The 
Hamiltonian is 



lit 



(1) 



where Vy is the hopping amplitude between sites i and 
j, a is the spin index, and the local interaction between 
electrons of opposite spin is U. We choose the hoppings 
such that the density of states becomes semi-elliptical, 
p(e) = v '4U 2 — e 2 / (27rV 2 ), and restrict our study to half- 
filling. Energy is measured in units of V and time in 
units of V -1 . To solve this model, we employ the dy- 
namical mean field theory (DMFT) ,±2 which has been 
extensively used to characterize the equilibrium phase 
diagram, and gives a qualitatively correct description for 
lattices with large coordination number^ The DMFT 
formalism has recently been reformulated for nonequi- 
librium situations ) 1 i 15 which gave insights into quench 
dynamics in the Hubbard model ) 16 ' and to nonequilib- 
rium states induced by external electric fieldsJ^— 

Nonequilibrium DMFT studies up to now have been 
restricted to the paramagnetic phases of the model, al- 
though the extension of DMFT to symmetry-broken 
phases is straightforward. 12 If the symmetry between 
spin-up and spin-down Green functions is not enforced, 
then the (single site) DMFT phasediagram for the half- 
filled, repulsive Hubbard model exhibits an antiferromag- 
netically ordered phase at low temperature (denoted by 
AFM in Fig. [J). For attractive U, one finds an anal- 
ogous phase diagram with AFM order replaced by s- 
wave superconductivity (at half-filling the superconduct- 
ing state is degenerate with a CDW phase, but in the 
doped system, superconductivity is more stable). The 
nature of the s-wave superconducting (or AFM insulat- 
ing) state changes qualitatively as \U\ crosses the value 
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FIG. 1: (Color online). Left panel: Antiferromagnetic phase 
diagram for the half-filled Hubbard model. The blue line with 
circles shows the exact DMFT result (from Ref. [2^), the red 
line with open circles the NCA approximation and the green 
line with stars the OCA approximation. Black diamonds in- 
dicate Teg(U) for interaction quenches from U = 4, T — 0.1. 
Right panel: NCA phase boundary as a function of filling at 
U = 8. The star indicates the effective temperature of the 
doped Hubbard model with the same magnetization as in the 
trapped state (see text). 



corresponding roughly to the maximum in the critical 
temperature. This is known as the "BCS-BEC" crossover 
in the cold-atom literature. Here, we will focus on 
the strongly interacting "BEC" (or localized moment / 
"Mott") regime, and study the stability of the ordered 
phase after a sudden switch of U into the paramagnetic 
regime. The quench provides an idealized but theoreti- 
cally well-controlled excitation procedure. Although it is 
artificial from the point of view of condensed matter ex- 
periments, we will see below that a qualitatively similar 
behavior of the long-time dynamics can be expected in 
photo-doping experiments. 

To study antiferromagnetic states within DMFT on a 
bipartite lattice, we have to solve impurity problems for 
each sublattice. For the semi-elliptic density of states 
(Bethe lattice), the hybridization function Aa, ct (Ab,<t) 
for the A (B) sublattice is given by the self-consistency 
condition A^ lCT = V 2 Gb,u (Ab, ct = V 2 Ga,<j), where 
G is the local lattice Green function. Together with 
the relation Aa,u = ^B,u (f° r pure Neel-type symme- 
try breaking), this leads to a single impurity calculation 
with self-consistency A CT = V 2 G S . To solve the nonequi- 
librium DMFT impurity problem, we use self-consistent 
strong-coupling perturbation theory. Figure [T] shows the 
equilibrium phase diagram obtained with the two lowest 
order implementations: the non-crossing approximation 
(NCA)2i and the one-crossing approximation (OCA). 25 
While NCA provides a rather good description of the ex- 
act QMC phase boundary 2 ^ in the strongly correlated 
regime, the maximum is shifted to f/^JSc ~ ^ (while 
L7max C w 4) and the phase boundary in the weakly cor- 
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FIG. 2: (Color online). Top panel: Time evolution of the 
magnetization for quenches from T = 0.1 and U = 4 to U = 6, 
7 and 8. The effective temperatures after these quenches are 
T eS = 0.113, 0.124 and T = 0.732, respectively, and the ar- 
rows indicate the corresponding thermal values of the magne- 
tization. Inset: same data on a different scale. Bottom panel: 
Time evolution of A(uj, t) for the minority spin after a quench 
from T = 0.1 and U — 4 to U — 8, and comparison to the 
thermal result (dashed). The pink curve corresponds to the 
spectral function of the long-lived nonthermal state. 



related regime is not well reproduced. OCA brings a sig- 
nificant improvement, with an almost correct position of 
the maximum and a quantitatively accurate description 
of the phase boundary in the strong correlation regime. 
Since we are mainly interested in quenches within the 
strongly correlated regime, we will use the NCA method 
for the real-time calculations. The techniques for the so- 
lution of the nonequilibrium DMFT equations and our 
implementation of the real-time NCA/OCA impurity 
solver have been explained in detail in Refs. H3 andl27l. 

We measure the time-dependent expectation values 
of the magnetization m — rif — n^, the double occu- 
pancy d = n^n^ (which yields the local energy E\ oc — 
Ud— fi(nf+ni)), and the kinetic energy per spin -Ekin.CT = 
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-iEijVijGf^t) (G<Jt,t>) = i(c{ a (t') C] At)))- 
The latter can be expressed (within DMFT) as a convolu- 
tion of the local Green function Gu^ a and the hybridiza- 
tion function Aj jCr . With Ej iCT denoting the self-energy 
for site i, we can write the lattice and impurity Dyson 
equations as [i9 t +/x-S i;0 .]Gy )CT -J2k v ikGkj,a 
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a Gii^ a , and with sublattice indices A and B one ob- 
tains 

£ kin , CT = -(i/2)[A At<7 G A ,a + A B , ff G B , CT ] < (t, t). (2) 



III. RESULTS 
A. Magnetization and spectral function 

The top panel of Fig. [5] plots the time-evolution of 
the magnetization for quenches from an initial equilib- 
rium state at U = 4 and T — 0.1, which is located deep 
inside the antiferromagnetic phase, to final states with 
U = 6, 7, 8. The sudden increase of U causes a change 
in the total energy E tot = Ei oc + ^kin.a of the sys- 
tem. E to t stays constant after the quench and allows us 
to define an effective temperature T c s(U) corresponding 
to the temperature of the thermal state with interaction 
U and total energy i?tot- This effective temperature is 
indicated by the black diamonds in Fig. [TJ While T e g 
initially decreases with increasing U due to the decreas- 
ing slope of constant entropy curves in this part of the 
phase diagram, it increases above U ss 4.5 and crosses 
the AFM phase boundary slightly above U = 7. After 
the quench to U = 8, the system is expected to thermal- 
ize in a paramagnetic (PM) state at high temperature 
(T ctf (U = 8) = 0.732). However, as shown in Fig. H af- 
ter a modest decrease the magnetization remains stuck 
for t > 20 at some large, nonthermal value, with no fur- 
ther relaxation evident on the timescales accessible in our 
simulations. Similarly, the quenches to U — 6 and 7 lead 
to a trapping in a nonthermal state, as is evident from 
the comparison to the thermal value of the magnetiza- 
tion (arrows). While the transient dynamics also shows 
interesting behavior, in particular slow oscillations in the 
amplitude of the order parameter (with period ~ 10), 
with superimposed rapid l//7-modulations (see inset of 
Fig. [2]), we will focus in this paper on the trapping phe- 
nomenon and the nature of the long-lived nonthermal 
state. 

To further characterize this state, we compute a time- 
resolved "spectral function" A(u,i) from the Fourier 
transform of the retarded Green function G let (t,t') = 
-ie{t-t>){{c(t),J{t>)}): 
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dt'e luj(t '- t} G rct (t',t). (3) 



This function is plotted for the minority spin in Fig. [21 
for the quench from U = 4, T = 0.1 to U = 8. While the 
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FIG. 3: (Color online). Time evolution of the double occu- 
pancy for the quench from U = 4, T = 0.1 to U = 8, 3, 2.5, 2 
(clockwise from the top left) and comparison to the thermal 
result (black arrow). The right panels show a narrow window 
around the thermal value and (as dashed line) an exponential 
fit to the long-time behavior. The left-pointing arrows in the 
left panels show the inital value of the double occupancy. 



red curve (t = 0) should not be confused with the spec- 
tral function of the initial equilibrium state, it exhibits 
Hubbard bands with pronounced spin-polaron peaks, as 
is typical for the AFM insulator^ Consistent with the 
left panel of Fig. [2l the spin imbalance shrinks from t = 
to t « 18 and then becomes time-independent for larger 
t. The spectral function of the trapped non-thermal state 
still features spin-polaron peaks, in contrast to the ther- 
mal spectral function at T e g = 0.732 (dashed lines). We 
conclude that despite the strong excitation of the system 
and the large amount of energy injected by the quench 
from U = 4 to U = 8, the antiferromagnetic order does 
not melt rapildy - instead, the system is trapped in a 
long-lived state with large magnetization, and with the 
typical spectral features of a magnetically ordered state. 



B. Double occupancy 

To understand the nature of this trapped state and the 
reason for its robustness, we consider the time-evolution 
of the average number of doubly occupied sites. The 
top left panel of Fig. [3] shows the result for the quench 
from U — 4, T — 0.1 to U = 8. There is a complicated 
transient regime (with 1/U oscillations) up to t w 20 
in which the double occupancy decreases from about 
0.044 to 0.021. At longer times, the double occupancy is 
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stuck at 0.021, even though the expected thermal value is 
0.0115 and thus almost a factor of two lower. This finding 
is consistent with theoretical and experimental studies of 
cold atom systems^ and with a previous investigation 
of pump-excited paramagnetic Mott insulators^ which 
show that the relaxation time of doublons grows expo- 
nentially with £/, leading - at large U - to a trapping in 
a nonthermal state characterized by a fixed number of 
doublons. 

Indeed, if one quenches to smaller interaction, the re- 
laxation time decreases and eventually becomes measur- 
able on the timescales accessible in the simulations (right 
hand panels of Fig. |3|). We have fitted the long-time 
behavior of the double occupancy with an exponential 
function decaying onto the thermal value and thus ex- 
tracted the relaxation times r w 400 for quenches from 
T = 0.1, U = 4 to U = 3.5, t w 47 (U = 3) and r « 7.8 
(U = 2.5). The [7-dendence of the relaxation time agrees 
with the analytical formula^ r = Ae a ^ u ^ los(u ^\ with 
A = 0.165 and a = 5.6. For the paramagnetic phase 
studied in Ref. [l9j, a smaller coefficient a was found, 
but due to the tails in the Gaussian DOS used there, 
the results cannot be directly compared by rescaling the 
hopping. 

Around U — 2 (close to the maximum critical tem- 
perature in NCA), the relaxation becomes so fast that 
the system already thermalizes within the time of the 
initial transient (bottom left panel of Fig. [3]). This find- 
ing of a qualitatively different relaxation pathway and 
fast thermalization at intermediate coupling is consistent 
with the results of Ref. [l^ for pump-excitations of the 
paramagnetic Mott insulator, and also with the result of 
Ref. [TH for interaction quenches from a noninteracting 
initial state. 



C. Comparison to a doped state 

The trapping of the double occupancy is well under- 
stood in the strong-coupling limit: For U 3> V , a unitary 
Schrieffer- Wolff transformation e s can be constructed 
order-by-order in V/U, such that successively terms of 
all orders in V/U are removed from the commutator 
[e~ s He s , d], and hence d = e s de~ s is conserved on ex- 
ponentially long times^ At second order, and after pro- 
jection to d = 0, one would obtain the t-J model, but 
here we encounter a more general situation, with d > 0. 
One can consider d as a rigorous definition of the number 
of free doublons, which differs from d by quantum fluctu- 
ations [S, d] + ... — 0(V/U) (even in the Mott insulator 
at T — 0, where d — 0). The initial very fast drop of 
d(t) on the timescale of a few inverse hoppings (Fig. [3]), 
is thus related to the reduction of quantum fluctuations 
due to the increase of U, or equivalently, a strengthen- 
ing of the local moments. Because the quench is faster 
than the timescale of a quantum fluctuation, quantum 
fluctuations are transformed into real doublon and hole 
excitations with a certain amplitude. This gives rise to a 




FIG. 4: (Color online). Comparison of the spectral function 
A(u,t — 18) of the trapped state after the quench from U = 
4, T = 0.1 to U — 8 (bold red line) to the majority-spin 
spectral function of a half-filled Hubbard model with T = 
0.0986 and to an infintesimally doped t-J model with J = 
4t 2 /U — 0.5 (top panel). The blue curve in the bottom panel 
corresponds to a doped Hubbard model with T = 0.0917 and 
chemical potential chosen such that the magnetization and 
double-occupancy of the trapped state are reproduced. (The 
t-J and doped Hubbard spectra are shifted by +4.2 and +2.4 
respectively.) 



nonzero d for t > 0, whose stability on exponentially long 
times prevents the system from thermalization. The ar- 
gument suggests that both the properties of the trapped 
magnetic states and the melting of the antiferromagnetic 
order can be related to the presence of injected free dou- 
blons and holes. To further support this fact it is useful 
to take a closer look at the spectral function. 

In Fig. [4] we compare the spectral function of the 
trapped state (t = 18) to equilibrium spectral functions 
of an infinitesimally doped t-J model with J = At 2 /U 
and identical magnetization, as well as to a doped AFM 
Mott insulator. The spectral function for one hole in the 
t-J model can be computed exactly within DMFTi 31 i 32 
The top panel of Fig. 2] reveals a good agreement of this 
t-J result and the spectral function of the half- filled Hub- 
bard model at [7 = 8, which shows that U — 8 is in the 
t-J limit, and that the NCA solution of the Hubbard 
model in this interaction regime is reliable. (At non- 
zero temperature, the spectral functions of the half-filled 
and infinitesimally doped Hubbard model are identical.) 
The temperature T = 0.0986 of the half-filled Hubbard 
model has been adjusted such that the magnetization 
corresponds to that of the trapped state, while the t-J 
spectrum has been rescaled by (1 + m)/2 = 0.815, and 
the t-J spectrum has been shifted on the frequency axis 
by +4.2. 

The figure also shows that the equilibrium results do 
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FIG. 5: Time-evolution of A < (uj,i) for a quench from U — 
4, T = 0.1 to U — 8. While the occupied states rapidly 
accumulate at the lower band edge in the AFM calculation, 
the occupation is nearly time independent in the PM case. 
The curves are slightly broadened by introducing a Gaussian 
factor e - "*-* - * ' in the integrand of Eq. @, which ensures a 
smooth cutoff of the upper integration limit. 



not very well reproduce the peak structures in the spec- 
tral function of the trapped state. The first peak near 
the gap edge is more prominent, while the higher en- 
ergy peaks are smeared out. According to the above 
strong-coupling argument, the trapped state might be 
close to an equilibrium state with additional free dou- 
blons and holes. Such a state with doublons and holes is 
not accessible in thermal equilibrium within the Hubbard 
model. In the paramagnetic case and at low energies, its 
properties arc different from a chemically doped state at 
the same doping^ Nevertheless, because the interaction 
of doublons with other doublons is very different from 
the interaction of doublons with holes, one may assume 
that the main correction to the spectral function at pos- 
itive energies (i.e., for inserting a doublon) is due to the 
presence of other free doublons, which motivates a com- 
parison of the spectral function of the trapped state to 
an electron-doped equilibrium state. We determined the 
doping level (n — 1.010) and temperature (T — 0.0917) 
of the equilibrium Hubbard model such that the mag- 
netization and double occupancy of the trapped state is 
reproduced. The corresponding point in the T-versus- 
filling phasediagram (right hand panel of Fig. [T]) is in- 
dicated by the black star symbol. The bottom panel 
of Fig. @] shows that the spectral function of this doped 
Hubbard model (shifted on the frequency axis by +2.4) 
is much closer to that of the trapped state. In particu- 
lar, the dominant first peak is well reproduced, while the 
damping of the peaks at higher energies is still not strong 
enough. This may be due to the missing interaction with 
holons. 



D. Occupation function 

Our results suggest that the trapped state is essentially 
an AFM Mott insulator with simultaneous electron and 
hole doping, comparable to a state that can be prepared 
by "photo-doping". The doping level is determined by 
the density of trapped doublons and holons and the ef- 
fective temperature is determined by the trapped magne- 
tization. In particular, the effective temperature of the 
trapped state after a quench from U = 4, T = 0.1 to 
U = 8 is apparently below the Neel temperature y^ r °P cd 
of this photo-doped state. This remarkable fact must 
be the result of "entropy cooling": the reduction in the 
magnetization from 0.85 to 0.63 releases a lot of entropy 
and this in turn keeps the trapped state below its Neel 
temperature. Evidence for the cooling effect is provided 
by the occupied density of states, 

A<(u,t) = Im / dt'e^ ~ l) (c f (t')c(t)) , (4) 

* Jt 

which corresponds to a time-resolved photoemission spec- 
trum for a quasi-steady state (see Fig. [5]). Since doublons 
and holes are inserted locally, they initially have a large 
kinetic energy, which leads to a broad distribution of the 
weight of A < (uj,t) in the upper Hubbard band. While 
the magnetization is decreasing, the weight becomes con- 
centrated at low frequencies (near the lower band edge), 
indicating that doublons and holes are being cooled by 
exchanging energy with the spin background. This is 
in contrast to a quench in the paramagnetic phase, af- 
ter which the occupation is not redistributed in time, 19 
because within DMFT the spin background of the para- 
magnetic state is completely disordered (no short range 
correlations) and hence does not allow energy exchange 
with doublons and holes. 

If the interpretation of entropy cooling is correct, one 
should be able to observe the transition of the trapped 
state from below Tjf pe to above Tjy° pe in the nonequi- 
librium dynamics, by changing the initial density of dou- 
blons and holes (T d r oped decreases with doping, see right 
hand panel of Fig. [T]). Indeed, as illustrated in Fig.[6l the 
magnetization decays to zero exponentially for quenches 
from [/i ni tiai < 3.3, T = 0.1 to U = 8, while the sys- 
tem is trapped in magnetized states for quenches from 
initial ^ 3.3, T = 0.1. A smaller initial U means a 
larger number of doublons and holes in the initial state, 
and thus also a larger doping of the effective photo-doped 
state. In other words, around [/initial = 3.3, the effective 
photo-doped state crosses Tjv- But even for [/initial < 3.3, 
where the AFM order decays to zero after the quench, the 
system is still trapped for very long times in a nonthermal 
state, as can be seen for example by comparing the double 
occupancy in the trapped state to the thermal value (top 
right panel of Fig. [B]). The trapped state still corresponds 
to a photo-doped state, but now at an effective temper- 
ature above TjJ° ped . In this context is also interesting to 
look at the timescale r for melting of the antiferromag- 
netic order, which diverges like r _1 cx U — U c close to 
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FIG. 6: Time-evolution of the magnetization (top left panel) 
and double occupancy (top right panel) for quenches from 
indicated initial U, T — 0.1 to U = 8. The bottom panel 
shows the inverse of the relaxation time r as a function of 
initial U. 



the transition (bottom panel of Fig. [6]) . This divergence 
resembles the critical slowing down at the equilibrium 
phase transition, which is related to the divergence of 
correlation times in equilibrium. 



E. Implications for superconducting states 

Finally, let us comment on implications of these results 
for other symmetry broken states, in particlular s-wave 
superconductivity (SC). As mentioned in the introduc- 
tion, the AFM phase of the half-filled repulsive Hubbard 
model can be mapped onto the s-wave SC phase in the 



half-filled attractive model. The mapping is a particle- 
hole transformation for the up-spin on the bipartite lat- 



tice: 



for i G A, 



Sit 



-c]^ for i e B, 



where A 



and B are sublattice indices. 34 Applying this transforma- 
tion to the repulsive model yields an attractive Hubbard 
model with U = —U. A photo-doped antifcrromagnet 
with doublons and holes transforms into a superconduct- 
ing or charge-ordered state with an enhanced number of 
unpaired electrons. Our finding of long-lived nonthermal 
AFM states in the repulsive model therefore implies the 
existence of similar long-lived nonthermal SC states af- 
ter quenches within the strongly-interacting regime on 
the attractive side. The similarity of the trapped AFM 
state to a photo-doped Mott insulator in turn suggests 
that a similar trapped SC state should appear after a 
photo-induced pair-breakup in the strongly interacting 
attractive regime. The stability of this nonthermal SC is 
then linked to the slow recombination of unpaired elec- 
trons. 

While it is clear that our simple model (with s-wave 
superconductivity) and the single-site DMFT formalism 
are in many respects not appropriate to describe the com- 
plicated physics of underdoped cuprates, one can spec- 
ulate that some of the physics discussed in this paper 
may be at play in the recently reported experiment by 
the Cavalleri groupj 1 ^ where a photo-induced transient 
SC state in a cuprate material was found to be stable 
for several tens of pico-seconds, and thus much longer 
than the light-pulse which stimulated apical oxygens and 
(by a yet unknown mechanism) induced the SC state. 
Long-lived, nonthermal symmetry broken states which 
are stabilized by the exponentially long life-time of dou- 
blons should also affect the relaxation dynamics in other 
strongly correlated, long-range ordered systems, such as 
CDW compounds. 
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